%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Study the temperature dependence of HMC by investigation fluctuations
% of function values (at the end of trajectories).
% 
% this function plot standard deviation of function value as a function of
% temperature
% 
% For example use parameters such as:
%
%   func_name = shekel
%   temperature range = as specified below
%   Nmd = 10
%   Nmc = 2500
%   NEachTempRun = 10
%   Pacc = 0.6
% 

% specify function name

func_name = 'shekel';

% read in HMC_param.used

[Tstart Tstep Tend Pacc Nmd Nmc Nruns] = read_param_used();


% Used that Nmc_to_equilibrium steps are needed before 'equilibrium' is
% reached

Nmc_to_equilibrium = 100;

if Nmc <= Nmc_to_equilibrium
   error('Nmc <= Nmc_to_equilibrium in code'); 
end


temperatures = Tstart:Tstep:Tend;


% initialise array containing s.d. of Epot values (at the end of
% trajectories)

sd_of_Epot_end = zeros([length(temperatures) Nruns]);  
sd_final = zeros([length(temperatures) 1]); 


index = 1;
for temp = Tstart:Tstep:Tend
  for i = 1:Nruns
    load_string = [func_name '/HMC_diagnosis' num2str(Nruns*(index-1)+i) '.out'];  
    
    fid = fopen(load_string, 'r'); % open file
    
    fgetl(fid); % read line and ignore it 
    fgetl(fid); % read line and ignore it
    
    x_sum = 0.0;  
    x_sq_sum = 0.0; 
    for i_Nmc = 1:Nmc
      fscanf(fid, '%f', 2); % read the first two number and ignore them
      if i_Nmc > Nmc_to_equilibrium
        x = fscanf(fid, '%f', 1);
        x_sum = x_sum + x;
        x_sq_sum = x_sq_sum + x^2;
      else
        fscanf(fid, '%f', 1);
      end
      fscanf(fid, '%f', 4); % read the remaining 4 numbers and ignore them
    end
    
    fclose(fid);
    
    N = Nmc - Nmc_to_equilibrium; % number of Epots used
    
    sd_of_Epot_end(index, i) = sqrt( x_sq_sum/N - (x_sum/N)^2 );
   
  end
  
  sd_final(index) = sum(sd_of_Epot_end(index, :))/Nruns;
  index = index + 1
end


plot(temperatures, sd_final,'x')
xlabel('T')
%xlabel('\langle E_{kin} \rangle')
%ylabel('sum(s.d. of \chi^2_{end})')